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S . Two  procedures  for  accelerating  the  convergence  of  the  SSOR  method 
are  considered,  one  based  on  conjugate  gradient  acceleration  and  the 
second  based  on  the  use  of  Chebyshev  acceleration.  Two  versions  of 
conjugate  gradient  acceleration  are  considered- -the  nonadaptive  and  the 
adaptive. 
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1.  Introduction.  In  this  paper  we  consider  the  symmetric  SOR  method 
(SSOR  method)  for  solving  the  linear  system 

(1.1)  Au  = b 

where  A is  a given  real  NxN  matrix  which  is  symmetric  and  positive  definite 
and  where  b is  a given  real  Nxl  column  matrix.  We  are  primarily  Interested 
in  cases  where  the  matrix  A is  large  and  sparse. 

The  basic  properties  of  the  SSOR  method  are  summarized  in  Section  2.  A 
procedure  is  given  for  estimating  the  optimum  relaxation  factor  a,  and  the  cor- 
responding spectral  radius  of  the  SSOR  matrix  S . These  estimates  are  in  tenns 

03 

of  the  spectral  radius  S(B)  of  the  Jacobi  method  and  another  quantity  ^ which 
can  usually  be  estimated  to  sufficient  accuracy. 

We  consider  two  procedures  for  accelerating  the  convergence  of  the  SSOR 
method.  The  first,  which  is  described  in  Section  3,  is  based  on  conjugate  gradient 
acceleration.  The  second,  which  is  described  in  Section  4,  is  based  on  the  use  of 
Chebyshev  acceleration.  To  carry  out  either  scheme  one  must  choose  the  relaxation 
factor  u'.  No  additional  parameters  are  required  to  apply  conjugate  gradient  ac- 
celeration. We  consider  two  versions  of  conjugate  gradient  acceleration:  the 
nonadaptive  version  where  the  value  of  a)  is  fixed;  and  the  adaptive  version  where 
one  chooses  a value  of  u)  and  improves  it  adaptively  based  on  the  results  of  the 
iterative  process. 
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In  order  to  apply  Chebyshev  acceleration  it  is  necessary  to  assume  values 

of  u and  S(fy  ).  We  consider  three  versions:  the  nonadaptive  version  where 
n 

neither  co  nor  the  estimate  of  S(»!J  ) varies;  the  partially  adaptive  version 
where  ai  is  kept  fixed  but  where  improved  values  of  determined 

adaptively;  and  the  fully  adaptive  version  where  improved  values  of  both  ai  and 
S(ty  ) are  determined  adaptively. 

Cl 

Preliminary  tests  have  been  carried  out  for  a very  simple  model  problem 
as  described  in  Section  5,  For  these  cases  the  adaptive  procedures  are  very  ef- 
fective in  the  sense  that  the  number  of  iterations  required  in  each  case  was  not 
substantially  greater  than  was  required  using  the  optimum  iteration  parameters. 

We  expect  that  the  accelerated  SSOR  method  with  adaptive  parameter  determina- 
tion will  prove  effective  in  a much  wider  class  of  problems.  Experiments  described 
by  Young  (2J  and  by  Benokraitis  [5]  indicate  that  the  nonadaptive  procedure  is 
effective.  Tests  on  similar  cases  based  on  the  adaptive  procedures  are  now  being 
carried  out. 

Acknowledgement . The  authors  would  like  to  acknowledge  the  contributions  of 
Mr.  Edward  Schleicher  in  preliminary  stages  of  the  work.  (See  [12].) 

Some  of  the  numerical  experiments  involving  the  Chebyshev  acceleration 
were  carried  out  by  James  Sullivan  in  the  Center  fcr  Numerical  Analysis  at  UT  Austin. 
We  would  also  like  to  acknowledge  the  help  of  Roger  Grimes  (L'T  Austin)  in  reading 
the  manuscript.  ~ 
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2.  The  basic  SSOR  method.  The  equation  (1.1)  can  be  written  in  the  form 


(2.1) 


u = Bu  + c 


where 


(2.2) 


/ 


B = I - D A 


c = D'^b. 


Here  D is  the  diagonal  matrix  with  the  same  diagonal  elements  as  A.  The  matrix 


B corresponds  to  the  Jacobi  method  of  iteration. 


The  matrix  A can  be  written  in  the  form 


(2.3) 


A = D - 


where  and  are  strictly  lower  and  strictly  upper  triangular  matrices. 


respectively.  Evidently 


(2.4) 


B = L + U 


whe  re 


(2.5) 


L . D-\  U . D-\ 


and  where  L and  U are  strictly  upper  triangular  and  strictly  lower  triangular 


matrices,  respectively. 
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The  SSOR  method  can  be  defined  in  terms  of  the  ordinary  (forward)  SOR  method 
and  the  backward  SOR  method.  The  ordinary  SOR  method  is  defined  by 


(2.6) 

where 


(n+1)  . (n)  , (F) 

u'  =ifu''+k 
a'  ai 


(2.7) 


£ ^ (I-a^L)  ‘^(a^U  + (1-0))!) 

a^ 


-1 


k^^^  = (I-coL)  V'C. 


a' 


The  backward  SOR  method  is  defined  by 


(2.8) 

where 


(n-tl)  (n)  , (B) 

u =‘?/u  +k 

U'  CO 


(2.9) 


= ( I - a^)  ■ ^ (a'L  I ( 1 -aO  I) 

u(B) 


(I  - a'U)  ^u'C . 
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liminate 
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we  obtain 

(2. 
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+1) 

S u 

05 

(n) 

4-  k 

O' 

.(n+1) 


where 


(2.12) 


A = V £ 

Cl'  a'  a' 

'^CD  ~ n'(2-u')  (I-u'U)  ^(I-a'L)  ^c, 


Mm sum 


i^iij  flip 


•"mil 


It  is  easy  to  show  that 


(2.13) 


dJ  = I - a.(2-a))  (I-ailO’^(l-tL'L)‘^D'^A 
a‘ 


I . liO;  ( i^D-  C )'^D(  -D-  C )‘^A 
a-.  u)  U (JD  L 


= I - Q A 


where 


(2.14) 

and  where 

(2.15) 


Q = cW^W 


O) 

2-a^ 


(2.16) 


1. 


Our  analysis  of  the  convergence  properties  of  the  SSOR  method  will  be  based  on 


chat  given  in  Young  [1974].  Let  m(B)  and  M(B)  be  the  smallest  and  largest 


eigenvalues  of  the  Jacobi  matrix  B.  We  remark  that  m(B)  < 0 and  M(B)  >0. 


Furthermore,  M(B)  < 1 since  A is  symmetric  and  positive  definite.  In  addition. 


the  eigenvalues  of  the  SSOR  matrix,  tS  , are  real,  non-negative,  and  less  than  unity. 


Let  P,  M,  and  m be  numbers  such  that 


m(B)  > m > -2  /P 

M(B)  < 2 ^ 

M < 1 

S(LU)  < P 


J 


(2.17) 


6 


It  can  be  shovm  that  S(B)  < 2 ; hence  if  we  have  a bound  M for  M(B)  such 

that  M > 2 we  can  replace  M by  2 , Similarly,  if  ra  <-2  /p’  we  can 

replace  m by  -2  . 


It  can  be  shown  that 


; 1 - ai(2-a') 


(2.18) 


I P > I or  if  P < 7 and  oj  < w* 

l-(iM-tu>  e ^ 4 - 


S(i  ) < <' 

O)  - I 


1 - a.(2-a') 


I-clto-Hd  P 


■2Z~  if  ^ 4 “ >0.'*, 


Here  for  P < T we  define  o)*  by 


(2.19) 


v/l-4p' 


Moreover,  the  bound  (2.18)  is  minimized  if  we  let 


(2.20) 


1 + \/i-ai+4p 


1 + n/1-4P 


if  M < 4p 


= ai*  if  M > 4p. 


The  corresponding  bound  on  S(i'  ) is  given  by 
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(2.21) 


S(.8  ) < < 

a^i  I 


'1-34  »4p 


, if  M ^ 4p 


1 - n/1  - 
1 t yi-4p 


= c*  - 1,  if  M > 4p, 


and  the  eigenvalues  A of  .b  will  be  in  the  range  0 A < S(t8  ) . 

‘^1 

We  shall  refer  to  the  value  of  given  by  (2.20)  as  a "good"  value  of 

a'.  Of  course,  is  not  necessarily  the  true  optimum  value  in  the  sense  of 

minimizing  S(t!ii). 

For  example,  if  the  matrix  A is  derived  from  a 5-point  difference 
equation  corresponding  to  a self-adjoint  elliptic  partial  differential  equation 


(2.22) 


L[u]  = (Au^)^  + 


and  if  we  rewrite  the  difference  equation  in  the  form  which  corresponds  to  the 


Jacobi  method 


(2,23)  u(x,y)  = e^(x,y)u(x  + h,y)  + 62(x,y)u(x,y  + h)  + B3(x,y)u(x  - h,y) 


where 


+ B^(x,y)u(x,y  - h)  + T(x,y) 


(2.24) 


6^(x,y) 


63(x.y)  = 


A(x  +-2.y) 
S(x,y) 


A(x  --j.y) 
S(x,y) 


B2(x,y) 


B^(x,y)  = 


C(x,y  + y) 

S(x,y) 


C(x.y-|) 

S(x,y) 


r 


s 


I 

1 

i 

1 


! 

1 


aiul 


S(x,y)  - A(x+^.y)  + A(x-^,y)  + C(x,y+^)  + CCx.y-y)  -h^F(x,y) 
i(x.y)  - -ir  . 


t hm  tlu'  t ommlas  for  the  SSOR  method  are  given  by 


/ 


(n+S)  , ■. 

u (x.y) 


(2  26) 


,(n+l) 


(x.y) 


a){  3jU^”\x  + h.y)  + 32*^  (x ,y  + h)  + - h.y) 

+ 3^u^"^'’\x,y  - h)  + T(x,y)}  + (1  - (x.y) 

u)t3jU^”  ^\x  + h,y)  + (x.y  + h)  + - h.y) 

+ 3^«^"'*’'^\x,y  - h)  + I (x.y)}  + (1  - u))u (x.y) 


\ 

If  one  uses  the  natural  ordering^.  P can  be  calculated  as 

(2.27)  p maxlpjfx.y) iP^(x-h.y)  . P2(x-h.y) ] . p^(x,y) [p^(x.y-h)  4 P2(x.y-h) 1 ) . 

An  upper  bound  M for  M(B)  can  be  obtained  from  (6.58).  page  1037  of  Young  and 
Gregory  [1973].^  The  relaxation  parameter,  ic.  is  calculated  using  (2.20). 


Application  to  the  Model  Problem  P 


We  now  consider  the  model  problem  defined  by 


-1 

0 ^ y 2;  1 with  u 0 on  the  boundary.  We  will  refer 
We  use  a mesh  size  of  h where  h ^ is  an  integer. 

*Thc  "natural  ordering"  for  points  in  a region  is  given  by  (x',y')  follows 
(x.y)  if  y'  > y or  else  y'  y and  x'  > x.  Thus  the  points  of  R^  are  treated  row 

by  row  from  left  to  right  starting  at  the  bottom  row. 

^This  bound  may  not  be  good  enough  in  some  cases.  In  that  case  an  adaptive 
procedure  may  be  used. 


(2,28) 


U 4 U 
XX  yy 


in  the  unit  square  0 ^ x 2;  1. 
to  this  as  "Model  Problem  P," 


I 
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Front  (2.2  7)  we  get 

(2.29)  P I . 

From  (6.58),  page  1037,  of  [1]  we  get 

(2.30)  M(B)  - cos  7rh. 

From  (2.20)  we  have 

(2.31)  a = " 

1 • y2(l-M(B)) 

which  is  very  close  to  the  optimum  co  , namely  , for  the  SOR  method. 
Moreover,  by  (2.21)  we  have 


1 . / izmi 

(2.32)  S(rf  ) V — ■■■  ■ • 

"1  1 , / iiMBi' 

For  Model  Problem  P,  it  can  be  shown  (see,  for  instance.  Young  [1974])  that 


(2.33) 

S(LU) 

1 2 TTh 

4 2 

Using 

1 2 TTh 

4 2 

for 

we  get 

(2.34) 


and 


(2.35) 


1 


_2 

+ V3  sin 


Sis  ) < 

a',  - 


1 - — 



1 . 

/3 


sin 


TTh 

2 


sin 


TTh 

2 


The  SSOR  method  requires  about  twice  as  much  work  per  iteration  as  the  SOR 

method  and  is  about  half  as  fast  as  the  SOR  method.  This,  at  first  sight,  would 

seem  to  preclude  the  use  of  SSOR.  However,  since  all  of  the  eigenvalues  of  the 

matrix  s are  real,  nonnegative,  and  less  than  unity  Young  [1974],  it  is  possible  to 
a 

accelerate  the  basic  SSOR  method  using  semi-iteration  or  variable  extrapolation 
as  shown  below.  The  resulting  method  is  faster  by  an  order-of-magnitude  than 
the  SOR  method.  In  oruer  for  this  gain  to  be  possible  it  is  sufficient  that 
S(LII)  - 7 be  of  the  same  order-of-magnitude,  in  some  sense, 

•4 


as  1 - M(B)  . 
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3 . Coniugrt-e  Gradient  Acceleration 


Let  us  now  apply  conjugate  gradient  acceleration  to  the  SSOR  method  defined 

by  We  note  that  Wcti  W ^ is  a symmetric  matrix,  where  W is  given  by 

a-’ 

(2.16) . Moreover,  as  stated  in  Section  2,  the  eigenvalues  of  are  real,  non- 

negative, and  less  than  unity. 

We  define  the  SSOR-CG  method,  i.e.,  the  conjugate  gradient  acceleration 

* ' 

procedure  based  on  the  SSOR  method,  by 


(3.1) 


= p ,(y  . u^")}  . (1-p  ,)u("-l> 

n • 1 n • i n 1 i 


where  is  the  pseudo- residual  vector  defined  by 


(3.2) 


and  the  y's  are  given  by 


(3.3) 


^n,l  - ('-'^n-l^ 


-1 


(3.4) 


(W5^"\we3  5(">) 


BQ 


Thus  RQ  is  the  Rayleigh  Quotient  of  the  vector  Wb^'^^  and  the  matrix  WGW  V 

n 1 

The  given  by 


/• 


(3.5) 


Pi  - 1 


Fl  .W£ 

L ■ (Wb<"*^>, 


(n)  w5(n))  ^ 

W6<"'^>)  ^n 


-1 


The  formulas  can  be  derived  from  Che  analysis  of  Concus,  Golub,  and 
O'Leary  [6).  See  also  Young  [11]. 
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The  L'onjiigato  gradient  acceleration  procedure  has  the  remarkable  property 


III'-  III'' 

where  u^*^^  - u (u  Is  the  exact  solution  of  (1.1))  corresponds  to 

L*v» 

conjugate  gradient  acceleration,  and  error  vector  corresponding 

to  any  semi-lterative  method  based  on  the  SSOR  method.  Note  that  for  any  vector  v 


we  have 


(3.8) 


(v,Av) 


Thus  Itvll  L could  he  computed  without  too  much  difficulty, 

a’" 

- 1 i 

Since  W(I-c>l  )W  (where  W is  given  by  (2,16))  and  A (!-,«  )A  are 

ui  in 

I 

symmetric  and  positive  definite,  we  could  use  either  W or  A for  W in 

, - lb  1 i 

(3.4)  and  (3.5).  The  choice  W D (— D - C ) minimizes  the  A -norm  of  the 

(0  U 

error,  and  the  choice  W = A minimizes  some  other  norm  of  tlte  error  at  each  step. 

An  inefficient  but  attractive  scheme  wovild  be  to  use  W D^(~D  - C,,)  in  the 

ii'  II 

iterative  and  adaptive  procedures  and  then  use  the  A -norm  for  the  stopping  pro- 
cedure. One  would  then  be  using  the  norm  whicli  is  being  minimized  in  the  stopping 
tests.  This  would  involve  an  extra  matrix/vector  opei'ation. 


I '3 


Kecurslve  Pscnido-  resldwals 


We  now  describe  an  alternative  formulation  of  the  SSOR-CG  method  defined 
by  (3.1) -(3.5).  This  alternative  formvilation  was  brought  to  our  attention  by 
Alan  Cline  of  The  University  of  Texas  at  Austin.  It  has  been  used  by  Dr.  Gene 
Golub  of  Stanford  University  and  others. 

In  the  formvilation  (3.1) -(3. 5)  one  is  required  to  do  two  matrix/vector 
multiplications  involving  S each  iteration,  namely,  u^*’^  in  (3.2)  and 
fS  in  (3.4).  We  can  reduce  this  to  one  such  matrix/vector  multiplication 

by  compviting  the  recursively.  Thus  by  (3.1)  and  (3.2)  we  have 


(3.8') 


ci(ti'l)  r o kC*!)  /)  /I 

n 1 1 n t L lo  nil  nil 


(n)  ^(n)  ^(n-1)  (n) 


rhe  calculations  could  be  done  as  follows.  Given  S b ,6  ,5 


to 


u^*'  and  p , we  could  compute  the  following: 
n 


Yn.l 

by 

(3.3) 

^^nil 

by 

(3.5) 

^^(n.l) 

by 

(3.1) 

,(ntl) 

by 

(3.8'). 

It  should  be  noted,  however,  that  this  procedure  involves  several  matrix/vector 
multiplications  involving  the  matrix  W.  These  can  be  eliminated  as  follows.  Let 

be  the  pseudo- residual  associated  with  the  SOR  method.  Thus  is  given  by 


(3.9) 


and  6 
(3,9') 


(n) 


^(n) 

is  given  by 

r(«) 


/ u 

U' 


(n) 


- u<") 


lO 


./  I u<">}  , - u<"^ 

a'  L\> 
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It  can  be 

shown 

that 

(3.10) 

W(I-rf  ) 

U' 

(3.11) 

Thus  we  have 

(3.12) 

(W5^"\' 

Moreover, 

the  A 

and 

(3.13) 

1 i 

7 d\i-/  ) 

g <1.' 

$ 


a' 

,(n) 


& 


,(n)  , .(n) 


] 


(n-1) 


5 = p, 


n+1  n+1 


(3-1^) 

The  p , and  y , are  given  by 
ntl  n+1  ° ^ 


(n) 


o)  a> 


" (l-Pn+l^^ 


(n-1) 


(3.15) 


j 1 , if  n = 0 


(3.16) 


' 'I  fi  . Isii  !_ 


-1 


, if  n > 1, 


A procedure  based  on  these  formulas  is  described  in  the  flow  chart  of 
Figure  3.1.  We  will  defer  discussion  of  the  stopping  tests  until  later. 


See  Hageman  and  Young  [l977],  Chapter  D. 


lb 


We  note  that  this  scheme  requires  the  application  of  the  matrix  £ 

U' 

to  5^  \ ‘^l  to  £ and  to  the  vectors  and  For 

Ci.' 

point  SSOR,  D is  a diagonal  matrix^  and  for  line  SSOR,  D is  tri-diagonal;  so 

this  requires  very  little  extra  work,*  Since  <tf  £ 6^"^  = iS  this  operation 

u'  u'  a' 

requires  as  much  work  as  one  full  SSOR  iteration.  Thus  the  total  work  per  iteration 
is  approximately  that  of  doing  one  SSOR  iteration. 

This  scheme  requires  considerably  less  work  per  iteration  than  that  given 
by  (3.1)  to  (3.5),  at  the  cost  of  increased  storage  requirements.  In  this  scheme, 
one  must  provide  additional  storage  for  the  A^"^  and  vectors. 


In  this  case  we  can  factor 
We  can  use  S 


T 

D into  the  form  S S where 
i 

instead  of  D throughout. 


S is  an  upper 


bi-diagonal  matrix. 
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Stopping  Tests 


Ideally,  we  would  like  to  require  that  our  final  approximation  u 
the  exact  solution  u of  (1,1)  satisfies 


(n) 


to 


(3.17) 


< c 


It  can  be  shown  (see  Hageman  and  Young  [l977j)  that  (3.17)  holds  if 


(3.18) 


V 2-a'  V 1-M 


(B) 


1 

1-S(dj  ) 


11^ 

11“  II 4 

D 


We  will  show  below  how  we  can  get  estimates  ^^(tS  ) of  S(«8^^^)  while  carrying  out 


the  SSOR-CG  method.  If  a good  estimate  M(B)  is  available  we  can  use 

the  stopping  test 


(3.18  ') 


where  is  the  SOR  pseudo- residual . This  follows  by  (3.11)  and  (2.15). 

In  the  adaptive  SSOR-CG  method  which  we  describe  later  in  this  section 
we  will  automatically  get  estimates  of  M(B) . If,  however,  we  are  not  using  the 
adaptive  procedure  and  if  we  do  not  have  a good  estimate  of  M(B),  then  we  can  be 
satisfied  with  requiring  that  our  final  approximatinn  u^*^^  satisfy 
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Since  our  iterative  process  yields  estimates  S («8  ) for  S(S  ) we  can  use  the 

E a'  a' 

test 

(3,20') 

u has  been  replaced  by  u^"^  in  (3.18')  and  (3.2![f)  . In  practice,  u is  not 
available,  and  as  n increases  ||u^  ^ ^ ^ become  a good  estimate 

of  ||u||  and  ||G||  . 

jj2  W 

We  remark  that  if  one  uses  either  (3.18’)  or  (3.20')  it  is  necessary  to 

compute  ||u^  ^ 11^’  ®fter  each  iteration.  However,  when  one  has 

reached  the  point  where  ||u^"^  - G||  ,,  or  ||u^")  - u|L , is  small,  it  is  not 

D ” 

necessary  to  change  ||u^"^  ||  or  ||u^"^  |L . In  practice,  ||u^"^  ||  and  ||u^"^  || 

D ” W 

are  not  changed  once  the  numerator  of  the  stopping  test  is  small. 
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(3.c«") 


For  our  estimate  S^(«8  ) we  use 
E u) 


S„(.S  ) - M(T  ) 

b u'  n 


where  M(T  ) Is  the  largest  eigenvalue  of  the  matrix  T . Here  T is  given  by 
n o o n n “ 


-d- 

Y 


-Yi)  y-(l'P2) 

r~  / 


1 

! 

1 V n V2 


(l-Y, 


(3.21)  T = 

n f 


/ 


/ 


(1-P3) 

Y2P2Y3P3 


(1-P3) 

Y2P2Y3P3 


-(1- 

Y, 


-Vl>  / -<^-^n> 
n-1  VVlVlVn 


L 


/ 


-(1-Pn) 


-(1-Y„) 


^n-l^n-l^n^n 


The  idea  of  using  the  eigenvalues  of  a tri-diagonal  matrix  derived  from  the 
conjugate  graduate  method  to  estimate  the  spectral  radius  of  an  iteration  matrix  is 
used  by  Concus,  Golub,  and  O'Leary.  See  also  Kaniel  [1966],  Paige  [1971],  and 
O'Leary  [1975]. 
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Adaptive  Determination  of  a'  for  the  SSOR-CG  Method 
So  far,  in  our  discussion  of  the  SSOR-CG  method  we  have  assumed  that  u'  is 
fixed.  Actually,  since  S(iS  ) is  a very  slowly  varying  function  of  a',  it  is  not 

U’ 

essential  that  the  exact  optimum  value  of  u'  be  used.  It  would  not  be  unreasonable 
to  simply  guess  at  a value  of  u-,  say  1.8,  and  use  the  non-adaptive  process  described 

above.  The  sensitivity  of  S(rt  ) to  u'  is  much  less  than  is  the  sensitivity  of 

S (/  ) . 

CO 

However,  it  is  possible  to  improve  a'  adaptively  using  the  matrices 

described  above.  Instead  of  working  with  a',  we  actually  work  with  which 

is  an  estimate  for  M(B) . We  assume  that  we  have  an  upper  bound  p for  S(LU) . 

(Such  a bound  should  not  be  much  greater  than  1/4  if  the  method  is  to  be  effective.) 

We  choose  an  estimate  Mj-(B)  such  that  M(B)--if  nothing  better  is  available 

let  M (B)  - 0.  We  then  compute  a>,  based  on  M (B)  and  p,  by  (2.20)  (with  ^ 

E E 

replaced  by  M^(B)>.  We  also  compute  Sj,(f8^)  by  (2.21). 

Before  beginning  with  the  adaptive  process  we  choose  an  "adaptive  factor" 

F such  that  0 < F < 1,  A typical  value  of  F is  3/4,  The  effectiveness  of 
the  process  is  quite  insensitive  to  F.  If  we  use  too  large  a va  ue  of  F,  say 


F very  close  to  unity,  we  will  change  parameters  very  often,  which  is  inefficient. 
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With  too  small  a value  of  F we  will  not  change  parameters  often  enough  and  we  will 
be  working  with  a poor  value  of  cn. 

After  each  iteration  we  compute  S'  - M(T^),  and  we  change  ^^.(8)  if 


(3.23) 


^1 

— < F 


where 


(3.24) 


= -log 


®(S„(.8  )) 

t a' 

fs  (s  )] 
0 — g- 
V S'  / 


(3.25) 


X2  = -log  '1>(S') 


Here  we  define  4’(x)  for  x e [0,1]  by 


(3.26) 


‘^(x)  = 


1 - /l^ 

1 + /T^ 


This  procedure  is  based  on  the  fact  that  Is  asymptotic  average  rate  of 

convergence  of  the  SSOR-SI  method  based  on  the  use  of  S„(e8  ),  if  the  true  spectral 

. ' E a' 

radius  S(t5^)  is  equal  to  S',  and  is  the  asymp.totic  average  rate  of  conver- 

gence  of  the  SSOR-SI  method  based  on  the  use^qf  ) = S',  'qgain  assuming  that 

the  true  spectral  radius  is  equal  to  S'  .•  ta^lo  of  these  convergence  rates 

is  at  least  F,  then  we  consider  that  the  convergence  is  fast  enough  and  we  do  not 
change  parameters.  If  the  ratio  is  less  than  F,  then  we  consider  that  the  con- 
vergence is  too  slow  and  we  change  the  parameters  to  and  S (^  ) in  order  to  obtain 

E u' 

faster  convergence. 

If  P < 1/4  and  if  at  any  stage  of  the  adaptive  process  we  have 


(3.27) 


-log  ^(a^*-!) 

-log  «>(S  (.8  )) 
E a'* 


then  we  can  settle  on  cd*  for  co  and  we  do  not  need  to  even  consider  changing 


r 


any  more.  If  (3.2  7')  is  satisfied  we  know  that  by  using  to  = to*  our  asymptotic 
average  rate  of  convergence  will  be  at  least  F times  that  which  we  could  estimate 
by  our  methods  for  any  other  to.  (It  should  be  noted,  however,  that  the  actual 
average  rate  of  convergence  for  some  to  might  be  greater  than  we  can  predict.) 

Having  decided  to  change  parameters  we  compute 


(3.28) 


^^^®^^new  = M’(B),  M^(B)) 


whe  re 


(3.29) 


«^(B) 


+ a 6)  - to(2-to) 


a)(ai  - (1+S')) 


(3.30) 


M^(B) 


, if  |(m(B))|  < M(B) 


^ , otherwise. 


New  values  of  u'  and  S (eS  ) are  obtained  by  (2.20)  and  (2.21),  respectively, 

£ CL) 

with  S . 

We  remark  that  the  value  of  M^(B)  is  obtained  by  setting  the  first  form 
of  the  right  member  of  (2.18)  equal  to  S'  and  solving  for  M.  We  note  that  the 
first  form  should  be  used  since  we  can  show  that  if  P < 1/4,  then  each  co  that 


we  use  will  not  be  greater  than  a^*. 


t t 

If  A has  Proparty  A or  is  an  L-matrix,  then  im(B)  < M(B) 
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Let  us  define  the  integer  s as  follows.  The  current  value  of  co  has  been 
(or  will  be)  used  for  the  first  time  in  the  computation  of  The  formulas 

for  the  adaptive  SSOR-CG  method  are  the  same  as  in  the  non-adaptive  case  except  that 

c 


(3.31) 


n+1 


V,  L 


1 ^ (d^a(">.d^a^">) 

■ Y„  .nK(n-l)  .^.(n-l) 


(D"A' 


,D"A' 


) 


-1 


if  n = s 


if  n > s 4 1 . 


The  overall  procedure  is  illustrated  in  Figure  3,2.  The  input  data  are 
u^^\.  M^(B),  F,  and  p.  We  set  n and  s equal  to  zero  and  compute  a^  and 

Next,  if 


S_(«S  ) as  well  as  the  initial  pseudo- residual  vectors  A^^^  and  6^^^ 
E u 


n ^ 0 we  find  M(T  ) where  T is  defined  by  (3.32).  (This  is  not  done  if  n = 0.) 
ri^  S T2 ^ s 

We  then  apply  the  stopping  test  (3.18')  with  S(<!S  ) replaced  by  S_(c8  ),  If  con- 

00  E cu 

vergence  has  not  occurred,  we  consider  whether  to  change  M (B) . 

£ 

The  "switch"  0 is  initially  set  equal  to  . However,  each  time  we  go 

to  we  test  whether  co*  would  be  satisfactory.  This  is  the  test  (3.27). 

If  at  any  time  (3.2  7)  is  satisfied,  then  we  set  to  , and  from  then  on 

M^(B)  and  co  are  not  changed. 

Assuming  that  0 

Mg(B)  should  be  changed.  If  Mg(B)  is  to  be  changed,  we  carry  out  the  change  using 

We  then  proceed 


we  then  use  the  test  (3.27)  to  determine  whether 


(3.28).  We  then  compute  new  values  of  a\,  Sg(t8^),  and  5^”^ 


to  calculate  jf  6^*^^  u P u and  We  are  then  ready 

05  n+1  n+1  ’ ’ 


for  the  next  iteration. 


The  computation  of  ;C  S 

to 


(n) 


’^n+l’  ''n+l’ 


^(n+1)^ 


, and  5 


(n+1) 


is  the 


same  as  in  the  nonadaptive  case  except  that  given  by  (3.31). 


The  matrix  T^  ^ is  used  instead  of  T^  since  when  we  change  ^(®)  (after 
the  s-th  iteration),  we  essentially  start  the  conjugate  gradient  process  over. 


hi 
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Accc’loratton  HuhihI  on  Sojiil -_1 


As  an  nlt«’rnatlvo  to  Lho  SSOK-Ctl  motliod  wo  can  vise  Ohobyshov  soml-1  torat  Ion 
to  accolorato  the  SSOK  motltod.  The  SSOK-Sl  method  Is  dotliu’d  hy  (.soo  Yoiiny,  |2h 


' I 1 I 1 t (1  - p ,>n' 

n I I n 1 1 11 1 I 


wlioro  th<’  psoiido-rostdnal  vi'ctoi:  is  given  hy  (1,1?)  and  wlu'te 


'll  I 1 ;?  - S (rii  ) 

ih 


IV,  (i-W‘ 


.1  I 

p , ( I - 70  p ) , 

11 1 1 4 n 


S(.H  ) 

O' 

In  order  to  apply  the  SSOK-SI  method  we  need  to  choosi'  .r  and  an  I’Stlmale 
S («H  ) for  S(<«  ).  Also,  i)  we  wish  to  use  the  stopping  test  (1.18')  w<>  neevi  an 

r.  ti'  O' 

estimate  Mj.(ll)  of  M(ll).  Dtlieivlse  we  use  file  stopping  test  (!1,.?0'). 

In  this  section  we  will  consider  two  adaptive  procedures.  The  llrsl,  which 
we  refer  to  ns  the  "partially  adaptive  procedure,"  Involves  fixing  o'  and  adaptively 
changing  the  estimate  S («s  ) for  S(tN  ).  The  other  procedure,  which  we  refer  to  as 

I'.  O'  <0 

the  "fully  adaptive  procedure,"  Involves  changing  o'  and  S (,«  ) adaptively.  I’he 

K O' 

procedure  where  o'  and  S (««  ) are  never  changed  Is  called  the  "non-adapt  1 ve  pro- 

I*,  tr 

cedure,"  I'or  the  fully  adaptive  procedure  we  work  with  M|.(ll)  and  hence  can 
use  the  stopping  test  (1.18'), 
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(4.11) 


1 - 


1 


1 


Here  F is  the  adaptive  factor  discussed  in  Section  3 for  the  adaptive  SSOR-CG 


procedure. 


Having  decided  to  change  S^(«S  ) we  let 

E u' 


(4.12) 


[S^(^  ))  1 
E a'  new 


max(S„(.8),  S'  (.8),  S;;(.8)) 

E u'  E a'  E O' 


where  Sl(t8  ) and  S"(t8  ) are  determined  as  follows.  To  get  S'{S  ) we  solve 

E U'  E a'  E o) 

the  Chebyshev  equation 


('•.13) 


115'=’  n, 


2r 


p/2 


1 + r^ 


1+r 


.•^iP 


where  r is  given  by  (4.11)  and 


(4.14) 


A 

r =• 


1 - A - is  (S  )/sl(s  )) 

E CL^  b g) 

1 + V^l-  (S  (.8  )/S'(.8  )) 

' E O)  E CO 


We  compute  S"  (i8  ) by  the  Rayleigh  quotient 
E a) 


(4.15) 


S''(.6  ) 
E ii> 


Having  computed  the  new  value  of  Sg(«8^)  we  let  s - n and  then  use 


(4,5) -(4.8)  to  continue  the  iteration  process. 
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The  overall  computational  procedure  is  illustrated  in  Figure  4,1.  As  in 
the  case  of  the  SSOR-CG  method,  we  determine  the  pseudo- residual  for  the 

SOR  method  for  each  iteration  as  well  as  the  pseudo- residual  6^*^^  for  the  SSOR 
method.  This  saves  matrix/vector  multiplications  by  W.  At  the  beginning  of 
each  iteration  we  compute  ^ and  . Usually  this  is  done  in  a straight- 

forward way,  but  if  n s+1  we  use  a special  procedure  as  described  below.  We 
then  carry  out  our  stopping  test.  We  use  (3,18')  if  a good  estimate  is 

available  for  M(B) . Otherwise  we  use  (3.20'). 

If  the  process  has  not  converged,  we  test  whether  S,(S  ) should  be 

L u' 

changed.  If  the  test  (4.9)  is  satisfied  or  if  n = 0,  we  compute  (S  («8  ))  = 

£ u'  new 

max(S„(rf  ),  S'(s  ),  S''(t8  ))  where  S’(S  ) is  the  solution  of  the  Chebyshev 
L u-  L CO  E Cl  E ui 

equation  and  S''(«8  ) is  the  Rayleigh  quotient  (4.15).  With  the  new  value  of 
E cr 

S_  («^  ) we  compute  p Y and  and  are  ready  for  the  next  iteration. 

E Cl  ntl  nii 
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We  use  a special  procedure  whereby  one  can  avoid  matrix/vector  multiplica- 
tions by  W and,  at  the  same  time,  can  compute  the  Rayleigh  quotient  S"(tS  ) each 

E o) 

time  changed  without  an  extra  matrix/vector  multiplication.  When 

we  decide  to  change  S„  (fB  ) we  compute  the  vector 

E to 


(^.16) 


~(n+l)  ^ ^(n)  ^ g(n)^ 


It  can  be  shown  that 


(^.17) 


= .8  + k - = .8  6<"> 


Thus,  for  S"(f8  ) we  have  by  (3.11)  (which  holds  for  and 

E to 


(^.18) 


s;'(.8  ) = rr i)rT 

E UN  (W&^  \w5^ 


(D^A<">,  D^A<‘'>) 


where 


(4.19) 


= £ ~(n+l)  ^ j^(F)  _ ~(n+l) 


We  note  that  to  get  the  Rayleigh  quotient  S"(«8  ) we  have  essentially 

E CO 

had  to  do  an  extra  SSOR  Iteration.  However,  we  can  recover  this  by  using  a special 

procedure  to  get  A^”^  and  6^*^^  on  the  next  iteration.  Thus  it  can  be  shown 
that  if  n = s+1  then 


^.ao) 


5^"^  = Y 5^"^  + (1-  Y 

n n 


A^"^  = Y + (1-  Y 
n n 


J 
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Fully  Adaptive  SSOR-Sl  Procedure 

So  far  we  have  been  considering  the  SSOR-SI  method  for  fixed  ai  where  we  make 

an  initial  estimate  ) for  S(t8  ) and  then  improve  S^(t8  ) adaptively. 

E a'  u)  E o) 

As  stated  in  Section  3,  S(^  ) is  a very  slowly  varying  function  of  cd,  and  it 

O) 

is  not  essential  that  the  value  of  ai  be  optimal.  Thus  it  would  not  be  unreason- 
able to  guess  at  a value  of  a',  say  1.8,  and  use  the  partially  adaptive  procedure 
described  above.  However,  as  we  show,  it  is  not  difficult  to  improve  a'  as 

well  as  S„(t8  ) adaptively. 

E CO 

When  we  let  a'  as  well  as  S„(«S  ) be  determined  adaptively  we  refer  to 

E ui 

the  procedure  as  "fully  adaptive.”  Actually,  as  in  the  case  of  the  adaptive 

SSOR-CG  method  considered  in  Section  3,  we  work  with  M (B) , an  estimate  of 

E 

M(B).  Also,  as  in  Section  3,  we  assume  that  we  have  an  upper  bound  3 for  S(LU) . 
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The  overall  procedure  for  the  fully  adaptive  SSOR-Sl  method  is  shown 

In  Figure  4.2.  To  begin  with,  we  estimate  a value  of  M (B)  such  that  0 < M (B)  < M(B) 

E E 

and  also  such  that  < 2\/^  . We  compute  cu  by  (2.20)  with  M = and 

S_(ei  ) by  (2.21). 

t U’ 

At  the  beginning  of  each  iteration  we  compute  and  5^”^ . We  then 

apply  the  stopping  test  (3.18').  If  convergence  has  not  occurred,  we  proceed  to  switch 


. Initially,  (^,  - . However,  if  at  any  time  a'*  can  be  shown  to  be 

satisfactory  (using  (3.27))  we  let  and  from  then  on  we  will  not  consider 

changing  M^(B).  We  then  apply  test  (4.9)  to  see  whether  ^^(B)  should  be  changed. 

If  so,  or  if  n = 0,  we  proceed  as  follows.  We  compute  a new  value  of  S (t8  ) as 

E a 

indicated.  The  value  of  S''(e8  ) given  is  the  Rayleigh  quotient,  as  was  shown  in  our 

t,  U' 

discussion  of  the  partially  adaptive  process.  Having  found  a new  value  of  S_(«S  ),  we 

E u'  ’ 

proceed  to  get  M_(B)  by  (3.28).  Then  we  compute  a new  ai  and  a new  S_(«8  ).  We 
t E u' 


then  compute  ^n+1’  ready  for  the  next  iteration. 
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5,  Numerical  results.  Several  experiments  were  run  for  the  SSOR  method 
using  the  conjugate-gradient  adaptive  scheme  and  using  the  semi-iterative  fully 
adaptive  scheme.  The  following  problem,  which  we  call  'taodel  problem  P, " was 


Numerical  experiments  on  more  general  problems,  such  as  those  considered  in  Young 
[l9  74]  are  now  being  carried  out. 

The  five-point  finite  difference  method  was  used  to  generate  the  matrix 

problem 

(5.2)  Au  = b 

for  mesh  sizes  of  h = 1/20,  1/40,  and  1/80.  For  these  cases  the  parameters  ^i, 

tDj^,  and  the  bound  for  ),  as  given  by  (2.30),  (2.31),  and  (2.32)  (with  P = 1/4), 

respectively,  are  shown  in  Table  1, 


TABLE  1 

VALUES  OF  OPTIMUM  PARAMETERS  FOR  MODEL  PROBLEM  P 


S(.8  ) 
0)1 


h = 1/20 

h - 1/40 

h = 1/80 

.98769 

.99692 

.99923 

1.  728  74 

1.85445 

1.92448 

.85451 

.52448 

.96151 

In  our  test  cases,  since  we  were  Interested  in  studying  the  effectiveness 
of  the  adaptive  schemes  rather  than  studying  effectiveness  of  stopping  procedures, 
we  generated  the  exact  solution,  u,  to  the  problem  (5.2)  and  iterated  until  the 
following  condition  was  satisfied: 


I 
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(5.3) 


K 


exact 


^ = 10-6. 

II  • II  ” ^ 


The  adaptive  procedures,  as  described  in  Sections  3 and  4,  were  used  for 
the  SSOR-CG  and  for  the  SSOR-SI  schemes,  respectively.  The  stopping  test  (5.3) 
rather  than  (2.18')  was  used  in  both  cases.  For  the  SSOR-SI  method  only  the 
fully  adaptive  scheme  was  used.  No  test  cases  were  run  using  the  partially  adap- 
tive procedure.  For  the  adaptive  process,  an  initial  guess  was  made  for 

p S(B).  It  is  necessary  that  1^,  and  in  the  absence  of  a better  guess 

£ 

can  be  set  equal  to  zero.  Initial  estimates  for  a'  and  S(t6  ) were  obtained 

from  (2.20)  and  (2.21)  using  the  estimated  parameter  . The  adaptive  procedures 

E 

of  Sections  3 and  4 were  then  used  to  correct  the  estimated  parameters. 

A non-adaptive  process  was  also  used  for  both  methods.  In  the  non-adaptlve 
process  a value  of  u'  was  chosen  and  not  changed.  In  the  case  of  the  non-adaptlve 

SSOR-SI  method  a value  of  S(«8  ) was  used  based  on  (2.21)  with  ^ replaced  by 

u.)  E 

“ “ 'if 

the  true  value  of  u and  with  p = 1/4, 

Table  2 gives  the  results  for  the  adaptive  and  non-adaptive  procedures 

where  a'  = .82  was  used  for  the  non-adaptive  process.  (This  corresponds  to  using 

(2.20)  with  1^;  = 0 and  p = 1/4.)  For  the  adaptive  cases  |a  was  initially  zero. 
L E 

The  value  F = 3/4  was  used  for  the  adaptive  cases. 

For  the  values  given  for  the  adaptive  SSOR-CG  method,  the  numbers  in  paren- 
theses correspond  to  a different  adaptive  scheme  where  at  least  two  iterations  were 
used  with  each  parameter  set.  The  numbers  look  somewhat  better  than  in  the  original 
adaptive  case. 


Clearly,  this  procedure  is  more  favorable  to  the  non-adaptive  SSOR-SI 
method  than  would  be  the  case  if  we  had  used  the  more  reasonable  procedure  of 
choosing  a value  of  p and  then  calculating  cd  and  a bound  for  S(i8  ) by  (2.18), 
letting  p = 1/4. 
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TABLE  2.  C(»1PARIS0N  OF  ADAPTIVE  AND  NON-ADAPTIVE  SCHEMES 


38 


Figure  5.1  shows  the  number  of  Iterations  required  to  satisfy  the  con- 
vergence criteria  (5.3)  for  the  non-adaptive  SSOR-CG  method  and  for  the  non- 
adaptive  SSOR-SI  method.  In  each  case  a value  of  cu  was  chosen  and  the  non- 
adaptive  procedure  was  used  as  described  above.  It  should  be  noted  that  the 
SSOR-CG  curve  is  much  flatter  than  the  SSOR-SI  curve,  and  also  lies  considerably 
below  the  SSOR-SI  curve.  This  indicates  that  the  SSOR-CG  procedure  is  much 
less  sensitive  to  choice  of  oi.  The  advantage  of  the  SSOR-CG  method  would  probably 
have  been  even  greater  if  we  had  not  used  the  true  value  of  n in  computing 
the  bound  on  S(t8  ) for  the  non-adaptive  SSOR-SI  scheme. 

U' 

The  results  presented  in  Tables  3 and  4 show  how  the  parameters  4 , a', 

E 

and  S(i9)„  changed  during  the  adaptive  process.  In  the  cases  considered,  the 
E 

initial  was  set  equal  to  zero.  It  is  evident  that  after  only  two  or  three 

L 

iterations,  good  values  of  the  parameters  are  available  through  the  adaptive 
process.  Furthermore,  even  the  iterations  occurring  before  the  selection  of 
the  final  parameter  set  are  not  "wasted"  iterations  since  the  iterant  u^*^^ 
is  being  improved,  although  not  as  much  as  would  be  the  case  with  the  optimum 
parameters . 
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tabu;  3.  SSOK-Si  FULLY 

> 

1 

< 

(F  3/4  in  all  cases) 

"‘f 

U' 

S(rf)  ) 
u'  E 

h 1/iX) 

Optimum  Value 

.98  769 

1 . 728  74 

.85451 

Iteration  No.  1 

0 

.828427 

.17157 

" " 4 

.9  7851 

1.6565  7 

.812156 

" " 15 

.98  742 

1 . 72621 

.85304  7 

" " 

Convergence 

h 1/40 

Optimum  Value 

.99692 

1.85445 

.92448 

Iteration  No.  1 

0 

.8284  2 7 

.17157 

" " 4 

.99034 

1 . 75597 

.8  7006 

" " 8 

.99633 

1.84218 

00 

" " clb 

Conve  rgence 

It 

Optimum  Value 

.99923 

1 , 92448 

.96151 

Iteration  No.  1 

0 

.828427 

.1715  7 

II  II  4 

. 99540 

1.82504 

.90852 

II  II  8 

.99896 

1.91280 

.95543 

II  II  j9 

ilonve  rgence 
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TABLE  4.  SSOR'CG  ADAPTIVE 
(F  = 3/4  in  all  cases) 


1 

1 

CO 

S(«8  )„ 

' CO  E 

i h = 1/20,  k = 1 

i 

j True  Value 

.98  769 

1.  728  74 

.85451 

Iteration  No,  1 

0 

.828427 

,1715  7 

i 

1 

M If  ^ 

.97851 

1.6565  7 

.812156 

? 

j 

If  If  ~j 

.98  742 

1. 72621 

.853047 

" " 16 

1 

Convergence 

i h = 1/40,  k = 2 

1 True  Value 

.99692 

1.85445 

.92448 

1 ! Iteration  No.  1 

1 

0 

.82847 

.17157 

! 

! ; " "4 

* 

,9913  7 

1.77713 

.8  71171 

" " 20 

Convergence 

h = 1/80,  k = 1 

1 

i 

True  Value 

.99923 

1.92448 

.96151 

‘ 

Iteration  No.  1 

0 

.82847 

.17157 

If  ff  2 

.99348 

1.  79  784 

.89047 

ff  If  7 

.99881 

1.91429 

.94801 

If  ff  22 

Convergence 
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Many  cases  were  run  for  various  values  of  ^ and  for  various  values 

of  F,  Figures  5,2  and  5.3  illustrate  the  effect  of  varying  p 5n  the  adap- 

E 

tive  procedure.  The  solid  lines  represent  the  non-adaptive  case;  i.e.,  is 

E 

chosen  and  the  parameters  are  calculated  based  on  and  are  not  changed 

E 

(u'  is  fixed  throughout  the  procedure)  . The  dashed  line  represents  the  fully 
adaptive  procedure;  here  is  specified  and  all  parameters  are  updated  when 

the  improved  parameter  set  is  calculated.  The  ordinate  represents  the  number 
of  iterations  required  for  convergence  and  the  abscissa  is  . The  graphs 

are  drawn  for  h = 1/40,  but  these  results  are  representative  of  all  of  the 
test  cases 

In  general,  lowering  the  value  of  F means  there  will  be  fewer 

parameter  changes,  and  setting  F closer  to  unity  has  the  effect  of  fine-tuning 
the  procedure.  In  this  case  there  will  be  more  parameter  changes  and  the  final 
set  should  be  closer  to  the  optimum  values.  However,  when  the  parameters  are 
changed,  the  convergence  rate  of  the  overall  method  drops  off  initially  and 
then  gradually  increases.  Thus,  changing  parameters  too  often  can  have  the  net 
effect  of  lowering  the  overall  convergence  rate.  In  choosing  F,  one  must  balance 
the  effect  of  changing  parameters  against  the  gain  which  will  be  realized  from 
improved  parameters.  Figure  5.4  is  a graph  of  the  number  of  iterations  required 
versus  the  parameter  F for  the  SSOR-CG  adaptive  procedure  and  for  the  SSOR-SI 
fully  adaptive  procedure.  In  general,  this  procedure  does  not  seem  to  be  too 


sensitive  to  values  of  F as  long  as  they  are  not  close  to  zero  or  to  one. 


SSOR-SI,  h = 1/40 


SSOR-SI  Fully  Adaptive 


Figure  5.4.  Number  of  Iterations  vs  F for  SSOR-CG  and  SSOR-Sl  Adaptive 
(h  = 1/40,  =0.,  p = 1/4) 
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6.  Summary  and  conclusions.  These  preliminary  test  cases  seem  to  indicate 
that  the  adaptive  procedures  presented  in  this  paper  are  effective  when  they  are 
coupled  with  both  the  SSOR-CG  and  the  SSOR-SI  methods.  It  has  been  shown  that  for 
this  problem  very  few  iterations  are  required  to  obtain  a good  set  of  parameters, 
and  these  preliminary  iterations  are  not  wasted,  in  the  sense  that  the  vector 
u^*'^  is  being  improved  on  each  iteration.  The  procedures  are  effective  even  when 
the  initial  guess  of  (i  is  the  worst  possible,"  i.e.,  = 0.  It  has  been 

shown  that  in  this  case  the  number  of  iterations  using  the  adaptive  scheme  is 
not  too  many  more  than  the  number  of  iterations  which  would  have  been  required 
if  the  optimum  parameters  had  been  known  from  the  start. 

These  results  are  preliminary  in  the  sense  that  numerical  experiments  have 
been  carried  out  only  for  the  model  problem  P.  We  are  now  testing  these  procedures 
on  a more  general  class  of  self-adjoint  elliptic  problems. 
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